home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / linalg / tridiag.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  6.8 KB  |  259 lines
  1. /* linalg/tridiag.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author: G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <stdlib.h>
  24. #include <gsl/gsl_errno.h>
  25. #include "tridiag.h"
  26. #include <gsl/gsl_linalg.h>
  27.  
  28. /* for description of method see [Engeln-Mullges + Uhlig, p. 92]
  29.  *
  30.  *     diag[0]  offdiag[0]             0   .....
  31.  *  offdiag[0]     diag[1]    offdiag[1]   .....
  32.  *           0  offdiag[1]       diag[2]
  33.  *           0           0    offdiag[2]   .....
  34.  */
  35. static
  36. int 
  37. solve_tridiag(
  38.   const double diag[], size_t d_stride,
  39.   const double offdiag[], size_t o_stride,
  40.   const double b[], size_t b_stride,
  41.   double x[], size_t x_stride,
  42.   size_t N)
  43. {
  44.   int status;
  45.   double *gamma = (double *) malloc (N * sizeof (double));
  46.   double *alpha = (double *) malloc (N * sizeof (double));
  47.   double *c = (double *) malloc (N * sizeof (double));
  48.   double *z = (double *) malloc (N * sizeof (double));
  49.  
  50.   if (gamma == 0 || alpha == 0 || c == 0 || z == 0)
  51.     {
  52.       status = GSL_ENOMEM;
  53.     }
  54.   else
  55.     {
  56.       size_t i, j;
  57.  
  58.       /* Cholesky decomposition
  59.          A = L.D.L^t
  60.          lower_diag(L) = gamma
  61.          diag(D) = alpha
  62.        */
  63.       alpha[0] = diag[0];
  64.       gamma[0] = offdiag[0] / alpha[0];
  65.  
  66.       for (i = 1; i < N - 1; i++)
  67.     {
  68.       alpha[i] = diag[d_stride * i] - offdiag[o_stride*(i - 1)] * gamma[i - 1];
  69.       gamma[i] = offdiag[o_stride * i] / alpha[i];
  70.     }
  71.  
  72.       if (N > 1) 
  73.         {
  74.           alpha[N - 1] = diag[d_stride * (N - 1)] - offdiag[o_stride*(N - 2)] * gamma[N - 2];
  75.         }
  76.  
  77.       /* update RHS */
  78.       z[0] = b[0];
  79.       for (i = 1; i < N; i++)
  80.     {
  81.       z[i] = b[b_stride * i] - gamma[i - 1] * z[i - 1];
  82.     }
  83.       for (i = 0; i < N; i++)
  84.     {
  85.       c[i] = z[i] / alpha[i];
  86.     }
  87.  
  88.       /* backsubstitution */
  89.       x[x_stride * (N - 1)] = c[N - 1];
  90.       if (N >= 2)
  91.     {
  92.       for (i = N - 2, j = 0; j <= N - 2; j++, i--)
  93.         {
  94.           x[x_stride * i] = c[i] - gamma[i] * x[x_stride * (i + 1)];
  95.         }
  96.     }
  97.  
  98.       status = GSL_SUCCESS;
  99.     }
  100.  
  101.   if (z != 0)
  102.     free (z);
  103.   if (c != 0)
  104.     free (c);
  105.   if (alpha != 0)
  106.     free (alpha);
  107.   if (gamma != 0)
  108.     free (gamma);
  109.  
  110.   return status;
  111. }
  112.  
  113.  
  114. /* for description of method see [Engeln-Mullges + Uhlig, p. 96]
  115.  *
  116.  *      diag[0]  offdiag[0]             0   .....  offdiag[N-1]
  117.  *   offdiag[0]     diag[1]    offdiag[1]   .....
  118.  *            0  offdiag[1]       diag[2]
  119.  *            0           0    offdiag[2]   .....
  120.  *          ...         ...
  121.  * offdiag[N-1]         ...
  122.  *
  123.  */
  124. static
  125. int 
  126. solve_cyc_tridiag(
  127.   const double diag[], size_t d_stride,
  128.   const double offdiag[], size_t o_stride,
  129.   const double b[], size_t b_stride,
  130.   double x[], size_t x_stride,
  131.   size_t N)
  132. {
  133.   int status;
  134.   double * delta = (double *) malloc (N * sizeof (double));
  135.   double * gamma = (double *) malloc (N * sizeof (double));
  136.   double * alpha = (double *) malloc (N * sizeof (double));
  137.   double * c = (double *) malloc (N * sizeof (double));
  138.   double * z = (double *) malloc (N * sizeof (double));
  139.  
  140.   if (delta == 0 || gamma == 0 || alpha == 0 || c == 0 || z == 0)
  141.     {
  142.       status = GSL_ENOMEM;
  143.     }
  144.   else
  145.     {
  146.       size_t i, j;
  147.       double sum = 0.0;
  148.  
  149.       /* factor */
  150.  
  151.       alpha[0] = diag[0];
  152.       gamma[0] = offdiag[0] / alpha[0];
  153.       delta[0] = offdiag[o_stride * (N-1)] / alpha[0];
  154.  
  155.       for (i = 1; i < N - 2; i++)
  156.     {
  157.       alpha[i] = diag[d_stride * i] - offdiag[o_stride * (i-1)] * gamma[i - 1];
  158.       gamma[i] = offdiag[o_stride * i] / alpha[i];
  159.       delta[i] = -delta[i - 1] * offdiag[o_stride * (i-1)] / alpha[i];
  160.     }
  161.       for (i = 0; i < N - 2; i++)
  162.     {
  163.       sum += alpha[i] * delta[i] * delta[i];
  164.     }
  165.       alpha[N - 2] = diag[d_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * gamma[N - 3];
  166.       gamma[N - 2] = (offdiag[o_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * delta[N - 3]) / alpha[N - 2];
  167.       alpha[N - 1] = diag[d_stride * (N - 1)] - sum - offdiag[o_stride * (N - 2)] * gamma[N - 2] * gamma[N - 2];
  168.  
  169.       /* update */
  170.       z[0] = b[0];
  171.       for (i = 1; i < N - 1; i++)
  172.     {
  173.       z[i] = b[b_stride * i] - z[i - 1] * gamma[i - 1];
  174.     }
  175.       sum = 0.0;
  176.       for (i = 0; i < N - 2; i++)
  177.     {
  178.       sum += delta[i] * z[i];
  179.     }
  180.       z[N - 1] = b[b_stride * (N - 1)] - sum - gamma[N - 2] * z[N - 2];
  181.       for (i = 0; i < N; i++)
  182.     {
  183.       c[i] = z[i] / alpha[i];
  184.     }
  185.  
  186.       /* backsubstitution */
  187.       x[x_stride * (N - 1)] = c[N - 1];
  188.       x[x_stride * (N - 2)] = c[N - 2] - gamma[N - 2] * x[x_stride * (N - 1)];
  189.       if (N >= 3)
  190.     {
  191.       for (i = N - 3, j = 0; j <= N - 3; j++, i--)
  192.         {
  193.           x[x_stride * i] = c[i] - gamma[i] * x[x_stride * (i + 1)] - delta[i] * x[x_stride * (N - 1)];
  194.         }
  195.     }
  196.  
  197.       status = GSL_SUCCESS;
  198.     }
  199.  
  200.   if (z != 0)
  201.     free (z);
  202.   if (c != 0)
  203.     free (c);
  204.   if (alpha != 0)
  205.     free (alpha);
  206.   if (gamma != 0)
  207.     free (gamma);
  208.   if (delta != 0)
  209.     free (delta);
  210.  
  211.   return status;
  212. }
  213.  
  214.  
  215. int
  216. gsl_linalg_solve_symm_tridiag(
  217.   const gsl_vector * diag,
  218.   const gsl_vector * offdiag,
  219.   const gsl_vector * rhs,
  220.   gsl_vector * solution)
  221. {
  222.   if(diag->size != rhs->size ||
  223.           (offdiag->size != rhs->size && offdiag->size != rhs->size-1) ||
  224.       (solution->size != rhs->size)
  225.           ) {
  226.     return GSL_EBADLEN;
  227.   }
  228.   else {
  229.     return solve_tridiag(diag->data, diag->stride,
  230.                          offdiag->data, offdiag->stride,
  231.              rhs->data, rhs->stride,
  232.                      solution->data, solution->stride,
  233.                      diag->size);
  234.   }
  235. }
  236.  
  237.  
  238. int
  239. gsl_linalg_solve_symm_cyc_tridiag(
  240.   const gsl_vector * diag,
  241.   const gsl_vector * offdiag,
  242.   const gsl_vector * rhs,
  243.   gsl_vector * solution)
  244. {
  245.   if(diag->size != rhs->size ||
  246.           offdiag->size != rhs->size ||
  247.           solution->size != rhs->size
  248.           ) {
  249.     return GSL_EBADLEN;
  250.   }
  251.   else {
  252.     return solve_cyc_tridiag(diag->data, diag->stride,
  253.                              offdiag->data, offdiag->stride,
  254.                  rhs->data, rhs->stride,
  255.                          solution->data, solution->stride,
  256.                          diag->size);
  257.   }
  258. }
  259.